20225421 HW5 Dahee Kwon

This homework is divided into 3 parts.
The first step is to convert light field images into 4D light field representation, and the second part is to refocus images at different depths. The last step is to compute an all-in-focus image and a depth map.
In this html file, I use the given chessboard light field image.

1.1 Initials

Load the light field image and create 5D array L(u; v; s; t; c).
u and v is the coordinates on the aperture
s and t is the coordinates on the lenslet array
c is the color channels.
Since given each lenslet covers 16 *16 pixels, the light field representation is 5D array with size (16,16,400,700,3).

lf_img = imread('./chessboard_lightfield.png');

imshow(lf_img);

lightfield = zeros(16,16,400,700,3);
 
for i = 1:16
    for j = 1:16
        for c = 1:3
            lightfield(i,j,:,:,c) = lf_img(i:16:6400,j:16:11200,c);
        end
    end
end

size(lightfield)
ans = 1×5
    16    16   400   700     3

1.2 Sub-aperture views

By rearranging the pixels in the light field image, we can create images with different camera aperture points (sub-aperture). This is equivalent to taking a slice of the light field 5D array with fixed u and v.
Since there is 16 number of u and v, we can make 16*16 sub-aperture images in total. Each of them is 400*700*3.

n = 0;

t = tiledlayout(16, 16, 'TileSpacing', 'none', 'Padding', 'compact');

for i = 1:16

for j = 1:16

n = n + 1;

nexttile;

img = uint8(squeeze(lightfield(i,j,:,:,:)));

imshow(img);

end

1.3 Refocusing and focal-stack generation

Combining parts of the light field can refocus the images at different depths.
If we average all of the sub-aperture views, then we can get the image with focused at a depth near the top of the chess board.

imgs = 0;

for i = 1:16

for j = 1:16

img = squeeze(lightfield(i,j,:,:,:));

imgs = imgs + img;

end

t = tiledlayout(1,1, 'TileSpacing', 'none', 'Padding', 'compact');

imshow(uint8(imgs/(16*16)))

If we want to move focusing points, then we should shift the sub-aperture images before averaging them.
To shift them properly, first we should shift the image (at fixed u and v) with respect to the center sub-aperture image.
The below code u and v is shifted version of u and v.

lensletSize = 16;
maxUV = (lensletSize-1)/2;
u = (1:lensletSize) -1 -maxUV;
v = (1:lensletSize) -1 -maxUV;

Then we can get re-focused images at different depths by shifting d*u (d*v) amount. d is the constant value that represents the depth.
Due to the inconsistency of the processing direction between lightfield and image (in matlab), note that we should make the shifting direction oppositely.
I use 8 number of depth points from {-1.50 to the 0.25}. (d=0 is equivalent to the non-shifted image)
If the front (bottom) of the chessboard is focused, the value of the depth should be smaller, so the d value representing depth appears negative in most cases.

S = 1:400;
T = 1:700;
d_list = {0.25, 0, -0.25, -0.5, -0.75, -1, -1.2, -1.5};
fstack = cell(8,1);
 
for d = 1:8
    imgs = 0;
    for i = 1:16
        for j = 1:16 
            img = squeeze(lightfield(i,j,:,:,:));
            
            s = round(S-d_list{d}*u(i)); s(s>400) = 400; s(s<1) = 1;
            t = round(T+d_list{d}*v(j)); t(t>700) = 700; t(t<1) = 1;
            img = img(s,t,:);
            imgs = imgs + img/(16*16);
        end
    end
    fstack{d} = uint8(imgs);
end

These are the refocusing images with each depth points.

t = tiledlayout(2,4, 'TileSpacing', 'none', 'Padding', 'compact');

for i = 1:8

nexttile;

imshow(fstack{i});

title(d_list{i});

end

1.3 All-focus image and depth from defocus

We can get a focal stack by merging the images into a new images.
In this process we should compute the pixel-wise weights and do weighted sum.
Weight is as follows

sigma1 = 2;
sigma2 = 0.5;
 
fstack_lum = cell(8,1);
fstack_low = cell(8,1);
fstack_high = cell(8,1);
wsharp = cell(8,1);

%1 Get luminance channel
for i = 1:8
    tmp = rgb2xyz(fstack{i});
    fstack_lum{i} = tmp(:,:,2);
end
%2 Create low frequency component
for i  = 1:8
    tmp = im2double(imgaussfilt(fstack_lum{i},sigma1));
    fstack_low{i} = tmp;
end
%3 Compute high frequency component
for i=1:8
    fstack_high{i} = fstack_lum{i} - fstack_low{i};
end
%4 Compute the (sharpness) weight
for i = 1:8
    tmp = im2double(imgaussfilt((fstack_high{i}).^2,sigma2));
    wsharp{i} = tmp;
end

wsharp is the weight per-pixels. To get wsharp, there are two parameters sigma1 and sigma2 which is the standard deviation of gaussian kernels.
I use various sigma1 and sigma2 then find the best one. (sigma1 = 2, sigma2 = 0.5)
As seen, although the effect of sigma is not significant (at least within the range I choose), if sigma values are increased, the sharpness of focused image is lost and the depth map also seems more blurry.
So I found parameters that focus on all pixels sharply enough, but the depth map is not too noisy. Note that I only conduct the experiments when sigma1>sigma2, because a sigma1 is used in the kernel for creating low frequency, and sigma2 is used to obtain the sharpness of the high frequency component.

Output

t = tiledlayout(1,1, 'TileSpacing', 'none', 'Padding', 'compact');

nom = 0; denom = 0;

for i = 1:8

nom = nom + wsharp{i}.*im2double(fstack{i});

denom = denom + wsharp{i};

end

focused_img = nom./denom;

imshow(focused_img)

When computing the depth map, negative d values are matter. Therefore, I add the offset (1.5) to make every depth value greater than 0 while keep the order.
It seems quite noisy but the depth information is included.

depthmap

nom = 0; denom = 0;

for i = 1:8

nom = nom + wsharp{i}.*(d_list{i}+1.5);

denom = denom + wsharp{i};

end

depthmap = nom./denom;

imshow(depthmap)

Bells & Whistles

To get blurred depth map, you should increase sigma values.
Also, normalizing weight values to handling extreme values is also one way to get better map.

nom = 0; denom = 0;

for i = 1:8

nom = nom + wsharp{i}.*(d_list{i}+1.5);

denom = denom + wsharp{i};

end

depthmap = nom./denom;

imshow((depthmap-min(depthmap,[],'all'))/(max(depthmap,[],'all')-min(depthmap,[],'all')));